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We present a calculation of the spectral functions and the associated rf response of ultracold fermionic atoms 
near a Feshbach resonance. The single particle spectra are peaked at energies that can be modeled by a modified 
BCS dispersion. However, even at very low temperatures their width is comparable to their energy, except 
for a small region around the dispersion minimum. The structure of the excitation spectrum of the unitary 
gas at infinite scattering length agrees with recent momentum-resolved rf spectra near the critical temperature. 
A detailed comparison is made with momentum integrated, locally resolved rf spectra of the unitary gas at 
arbitrary temperatures and shows very good agreement between theory and experiment. The pair size defined 
from the width of these spectra is found to coincide with that obtained from the leading gradient corrections to 
the effective field theory of the superfluid. 
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I. INTRODUCTION 

The existence of well defined, non-interacting quasiparti- 
cles above a possibly strongly correlated ground state is a 
central paradigm of many-body physics. In interacting Fermi 
systems, this concept applies both in a Fermi liquid and in a 
BCS-like superfluid state, whose elementary excitations have 
an infinite lifetime at the Fermi surface. More generally, the 
nature of quasiparticle excitations may be used to character- 
ize many-body ground states both with or without long range 
order [1]. Typically, it is only near a quantum phase transition 
between ground states with different types of order, where a 
quasiparticle description fails and is replaced by a continuum 
of gapless excitations [2]. In our present work, we discuss ul- 
tracold fermionic atoms with a tunable attractive interaction. 
The ground state is a neutral s-wave superfluid at arbitrary 
coupling. Thus, it has gapless bosonic quasiparticles of the 
Bogoliubov- Anderson type with a linear spectrum to = c s q. 
Its fermionic excitations have a finite gap. Within a BCS de- 
scription, the associated Bogoliubov quasiparticles are exact 
eigenstates of the interacting system at arbitrary momenta. As 
will be shown below, this central feature of the BCS picture 
of fermionic superfluids fails for the strong coupling situation 
that is relevant in the cold gases context, where the excita- 
tion energy is no longer exponentially small compared with 
the Fermi energy. In this regime, the fermionic particle exci- 
tations acquire a significant lifetime broadening even at zero 
temperature, except near the dispersion minimum (or maxi- 
mum for holes), where there is no available phase space for 
decay. The lifetime broadening arises both from the resid- 
ual interaction between quasiparticles and their coupling to 
the collective sound mode. Moreover, the particle-hole sym- 
metry characteristic for the Bogoliubov quasiparticles of the 
BCS theory is violated in the strong coupling regime. With 
increasing temperatures, the particle- and hole-like branches 
merge into a single broad excitation branch with a free parti- 
cle like dispersion, shifted by the binding energy. 

Fermions with a tunable attractive interaction and the as- 
sociated BCS-BEC crossover have been studied experimen- 
tally using ultracold Fermi gases near a Feshbach resonance 



[3, 4, 5]. The fact that the balanced system with an equal num- 
ber of particles in the two different hyperfine states ('spins') 
that undergo pairing is superfluid at sufficiently low temper- 
atures has been inferred from the observation of a finite con- 
densate fraction on the BCS side [6] and from the collective 
mode frequencies in a trap that agree with superfluid hydro- 
dynamics [7, 8]. It was demonstrated quite directly by the 
observation of a vortex lattice in the rotating gas, that evolves 
continuously from the BEC to the BCS side of the transition 
[9]. To study the excitation spectrum, in particular the evo- 
lution of the expected gap for fermionic excitations due to 
pairing, rf spectroscopy was performed by Chin et al. [10]. 
The interpretation of these measurements [11] in terms of a 
an effective 'pairing gap', however, is made difficult by the 
existence of strong final state interactions and the fact that the 
signal is an average over the whole cloud, with a spatially de- 
pendent excitation gap. For a homogeneous system, the aver- 
age rf shift is in fact dominated by large mean-field effects and 
final state interactions [12, 13, 14, 15] and is hardly changed, 
even if superfluidity is suppressed by a rather strong imbal- 
ance [16]. The problems associated with final state interac- 
tions and the inhomogeneity of the cloud have been overcome 
only recently by the possibility to perform spatially resolved 
rf measurements [17], combined with a suitable choice of the 
hyperfine states which undergo pairing and the final state of 
the rf transition [18]. Moreover, it has also become possi- 
ble to measure rf spectra in a momentum resolved way [19]. 
This opens the possibility to infer the full spectral functions, 
as suggested theoretically by Dao et al. [20]. 

Our aim in the following is to present a calculation of the 
spectral functions and the associated rf response of strongly 
interacting fermions which covers the whole regime of tem- 
peratures both above and below the superfluid transition and 
also arbitrary coupling constants. The theory is based on a 
conserving, so-called ^-derivable approach to the many-body 
problem due to Luttinger and Ward, in which the exact one- 
particle Green functions serve as an infinite set of variational 
parameters. This approach has been used previously to de- 
scribe the thermodynamic properties of the uniform [21] and 
the trapped gas [22]. The Luttinger- Ward formulation of the 
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many-body problem relies on expressing the thermodynamic 
potential ft[G] in terms of the exact Green function G. The 
condition that the functional 0[G] is stationary with respect to 
small variations of the Green functions then leads to a set of 
integral equations for the matrix Green function G which have 
to be solved in a self-consistent manner. Since the Green func- 
tions contain information about the full dynamical behavior 
via the imaginary time dependence of the Matsubara formal- 
ism, the Luttinger-Ward approach not only provides results 
for the equilibrium thermodynamic quantities but also deter- 
mines the full spectral functions upon analytic continuation 
from Matsubara to real frequencies. This is done explicitly in 
our present work, using the maximum-entropy technique. 

The paper is organized as follows: in Sec. II we introduce 
the Luttinger-Ward formalism and discuss the calculation of 
the momentum and frequency dependent spectral functions. 
The relation between the spectral functions and the experi- 
mentally measured rf spectra is outlined in general and dis- 
cussed in the BCS and BEC limit, where analytical results are 
available. We also discuss the behavior of the rf spectra at 
high frequencies and the associated contact coefficient intro- 
duced by Tan [23] and by Braaten and Platter [24]. In Sec. 
Ill, we show that a pair size can be defined via the momentum 
dependence of the superfluid response, in analogy to the non- 
local penetration depth in superconductors. Using an effective 
field theory due to Son and Wingate [25], we find that the re- 
sulting pair size of the unitary gas coincides with that inferred 
experimentally from the width of the rf spectrum [18]. The 
numerical results and the physical interpretation of spectral 
functions and rf spectra obtained within the Luttinger-Ward 
approach are discussed in Sec. IV, both in the normal and su- 
perfluid phase. These results are compared quantitatively with 
measured data. A summary and discussion is given in Sec. 
V. There are two appendices, one on the maximum-entropy 
method and one on a perturbative calculation of the quasipar- 
ticle lifetime due to interactions with the collective mode. 



a delta function in three dimensions leads to no scattering at 
all, the bare coupling strength 



00(A) = — 



9 



2aA/7r 



(2.2) 



needs to be expressed in terms of renormalized scattering 
amplitude g = Anh 2 a/m that is proportional to the s-wave 
scattering length a and an ultraviolet momentum cutoff A 
that is taken to infinity at fixed g. The limiting process 
<7o(A — > oo ) — > —0 accounts for the replacement of the bare 
delta potential by a pseudopotential with the proper scattering 
length. The description of a Feshbach resonance by a sin- 
gle channel Hamiltonian of the form given in (2.1) is valid 
for the experimentally relevant case of broad Feshbach reso- 
nances, where the effective range r* of the resonant interac- 
tion is much smaller than the Fermi wavelength [4]. 

We consider a homogeneous situation described by a grand 
canonical distribution at fixed temperature and chemical po- 
tential. The grand partition function 



Z = Tr{exp(-/3[if - fiN])} 
then determines the grand potential 

fl = Q(T,(jl) = -P^lnZ . 



(2.3) 



(2.4) 



For a quantitative discussion of the results, it is more conve- 
nient to switch to a canonical description at a given density n 
by a Legendre transformation to the free energy F = Q + fiN. 
Within our zero range interaction model, the Fermi system 
at total density n = kp/3ir 2 is then completely character- 
ized by two parameters: the dimensionless temperature 9 = 
ksT/eF and the dimensionless inverse interaction strength 
v = 1/kpa. In the special case of an infinite scattering length 
(the so-called unitarity limit), the parameter v drops out. The 
resulting spectral functions A(k, e) are then universal func- 
tions of 9 and the dimensionless momentum and energy scales 
k/kp and e/ef that are set by the density of the gas. 



II. LUTTINGER-WARD THEORY, SPECTRAL 
FUNCTIONS AND RF RESPONSE 



Our calculation of the spectral functions for a dilute system 
of ultracold fermionic atoms is based on a Luttinger-Ward ap- 
proach to the BCS-BEC crossover, that has been presented in 
detail previously [21, 26]. As a starting point, we use the stan- 
dard single-channel Hamiltonian, that contains the essential 
physics of the BCS-BEC crossover in a dilute gas of ultracold 
fermionic atoms with a short range (s-wave) interaction [4] 



H=J d 3 r^|^VV,+ (r)][V^(r)] 



(2.1) 



Here VvM and (r) are the usual fermion field operators. 
The formal spin index a labels two different hyperfine states 
which interact via a zero-range delta potential go S(r). Since 



A. Luttinger-Ward formalism 

In thermal equilibrium at temperature T the properties of an 
interacting fermion system which exhibits a superfluid tran- 
sition are described by two Matsubara Green functions, the 
normal Green function (T denotes the standard time ordering) 

(T[Vv(r, r )V+, (r', r')]) = S aa> G(v - r', r - r') (2.5) 

and the anomalous Green function 

(T[W(r, r)W (r', r')]) = e aa , T{v -r',r- r') (2.6) 

where the antisymmetric Levi-Civita tensor E aa i represents 
the spin structure of s-wave pairing. In the translation invari- 
ant and stationary case studied here, it is convenient to switch 
to a Fourier representation of the Matsubara Green functions. 
The normal and anomalous functions (2.5) and (2.6) can then 
be combined into a matrix Green function 
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with momentum variable k and fermionic Matsubara frequen- 
cies u) n = 2n(n + 1/2)/ (3h with n € Z. The nondiagonal 
elements represent the order parameter of the superfiuid tran- 
sition. Using the matrix Green function (2.7), it is possible 
to generalize the Luttinger-Ward formalism [27] to superfiuid 
systems [21, 26]. In particular, the grand thermodynamic po- 
tential (2.4) can be expressed as a unique functional of the 
Green function (2.7) in the form 

Q[G] =/3- 1 (-|Tr{-lnG+[(?o 1 G-l]}-$[G]) . (2.8) 

The interaction between the fermions is described by the func- 
tional $[G], which can be expressed in terms of a perturba- 
tion series of irreducible Feynman diagrams. The full matrix 
Green function G is then determined uniquely by the condi- 
tion that the grand potential functional (2.8) is stationary with 
respect to variations of G, i.e. 

m[G]/8G = . (2.9) 

It is important to note, that the thermodynamic potential Sl[G] 
is a functional of the exact Green function G. The formalism 



of Luttinger and Ward thus leads via (2.9) to a self-consistent 
theory for the matrix Green function. 

Since the exact form of $[G] is unknown, we employ a lad- 
der approximation [21, 26, 28]. In the weak coupling limit, 
this is exactly equivalent to the standard BCS description of 
fermionic superfluids. In the BEC limit, where the fermions 
form a Bose gas of strongly bound pairs, the ladder approxi- 
mation correctly accounts for the formation of pairs (i.e. the 
two-particle problem). The residual interaction between the 
pairs, however, is described only in an approximate manner. 
Indeed, it turns out [21, 28] that in the BEC limit the ladder ap- 
proximation for the functional $[G] gives rise to a theory for 
a dilute Bose gas with repulsive interactions that are described 
by a dimer-dimer scattering length add = 2a. This is a qual- 
itatively correct description of the BEC limit of the crossover 
problem, however from an exact solution of the four-particle 
problem in this limit the true dimer-dimer scattering length 
should be add = 0.6a [29]. 

The ladder approximation leads to the following closed set 
of equations for the matrix of single particle Green functions 
(2.7) [21, 26, 28] 



G~U^n) 
S aQ ,/ (k, u> n ) 



G 0,L'( k ' W «) _ £ <««'( k . w n) . 

^ G Q / Q (k + K, uj n + n n ) T aa , (K, tt, 



(2tt) 3 p 
d 3 k rl 



r-i,(K,0„) = &J y+J ^[^G aa ,(K--k,Sl n -w n )G aa ,(k,u n ) 



h 2 k 2 



(2.10) 
(2.11) 

(2.12) 



r 



Here, 

^O^a'^^n) - I q 







- [ivn + e k - n\ 

(2.13) 

is the inverse free Green function where = h 2 k 2 /2m. Fur- 
thermore, 



Si 



A 

A* 



(2.14) 



is a k- and uj n -independent matrix, whose off-diagonal ele- 
ments represent the order parameter of the superfiuid transi- 
tion. By definition, A is related to the anomalous Green func- 
tion T(k, r) by the renormalized gap equation 



A 
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d 3 k 



(2tt)= 



F(k,r = 0) + A 



h 2 k 2 



(2.15) 



The vertex function r QQ /(K, O n ) defined in (2.12) may be 
identified with the T matrix for the scattering of two particles 
in a many -body Fermi system. Since G aa > (k, ui n ) is the exact 
one-particle Green function, the vertex function is that of a 
self-consistent T matrix approximation. The Luttinger-Ward 



approach in ladder approximation is thus equivalent to a self- 
consistent T-matrix approximation. The specific structure of 
the GG term in (2. 1 2) with respect to the Nambu indices a and 
a' implies that the particle-particle ladder is considered here, 
which properly describes the formation of Fermion pairs in 
normal and superfiuid Fermi systems. 

As a result of the Goldstone theorem, a neutral superfiuid 
Fermi system must exhibit a gapless Bogoliubov-Anderson 
mode. Formally, this is guaranteed by a Ward identity, which 
can be derived from the Luttinger-Ward formalism for any 
gauge invariant functional $[G]. This functional defines an 
associated inverse vertex function which in short-hand nota- 
tion is given by 



r- 1 = rr x 



Y 



(2.16) 



where T\ = — <5 2 $[G]/<5G 2 is the irreducible vertex and 
X = —GG is the pair propagator. The existence of a Bo- 
goliubov Anderson mode is then guaranteed by the property 
that r _1 has an eigenvalue A(K, f2„) which has to vanish for 
K = and 0„ = [26]. This Ward identity is equiva- 
lent, in the present case, to the well known Thouless crite- 
rion [30]. Unfortunately, the inverse vertex (2.12) obtained 
from our self-consistent ladder approximation does not agree 
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with the exact inverse vertex function as defined by Eq. (2.16). 
As shown in our previous publication [21], however, the re- 
quirement of a gapless Bogoliubov-Anderson mode can be 
imposed on (2.12) as an additional constraint by choosing a 
modified coupling constant in the renormalized gap equation 
(2.15). This modified approach is still compatible with the 
Luttinger-Ward formalism so that our method is both conserv- 
ing and gapless. In the following numerical calculations we 
always employ this modified approach which is described in 
detail in Ref. [21]. 

In a homogeneous gas, the normal to superfluid transition 
is a continuous phase transition of the 3D XY type along the 
complete BCS to BEC crossover. By contrast, our approach 
[21] gives rise to a weak first-order superfluid transition be- 
cause the superfluid phase of the Luttinger-Ward theory does 
not smoothly connect with the normal-fluid phase at a single 
critical temperature 8 C = ksTc/sF- Fortunately, this prob- 
lem is confined to a rather narrow regime of temperatures. In 
particular, at unitarity, the upper and lower values for 9 C are 
0.1604 and 01506, which is within the present numerical un- 
certainties in the determination of the critical temperature of 
the unitary gas [31, 32]. For our discussion of spectral func- 
tions in the present work, which does not focus on the critical 
behavior near T c , the problem with the weak first order nature 
of the transition is therefore not relevant. 

Keeping these caveats in mind, the ladder approximation 
for the Luttinger-Ward functional provides quantitatively re- 
liable results for the thermodynamic properties of the BCS- 
BEC crossover problem [21, 22]. This applies, in particu- 
lar, for the most interesting regime near unitarity, where weak 
coupling approximations fail. As an example, the value of the 
critical temperature T c /Tp = 0.16 right at unitarity agrees 
with recent quantum Monte-Carlo results T c jTp = 0.152(7) 
for this problem within the error bars [31, 32]. It is also 
consistent with recent calculations of the onset temperature 
of a finite condensate density [33]. Moreover, there is also 
quite good agreement with field-theoretic results for ground 
state properties, which are characterized by a single univer- 
sal constant, the so called Bertsch parameter £(0) defined e.g. 
by ^{T = 0) = £(0)£z? at unitarity [4]. In fact, the value 
£(0) = 0.36 obtained within the Luttinger-Ward approach 
[21] agrees perfectly with the result £(0) = 0.367(9) from 
an e = 4 — d expansion up to three loops [34] and - in particu- 
lar - with the more recent value £ = 0.36 ± 0.002 obtained by 
Nishida [35]. Variational Monte Carlo calculations [36, 37] 
or a Gaussian fluctuation expansion around the BCS mean- 
field results [38, 39], in turn, give somewhat higher values 
£(0) = 0.42(1) or£(0) = 0.40, respectively. 



B. Spectral functions 

The Matsubara Green function G(k, u) n ) can be expressed 
in terms of a spectral function A(k, e) by using the Lehmann 
spectral representation [40] 



Q(k,w n ) = J de 



A(k,s) 



The spectral function associated with the normal, single par- 
ticle Green function Q(k,uj n ) is positive A(k,e) > and 
normalized according to 

J deA{k,e) = 1 . (2.18) 



It can be decomposed into two contributions 

A(k,e) = A + (k,e)+A-{k,e) 



(2.19) 



which describe the particle and hole excitation part of the 
complete excitation spectrum. The individual contributions 

A+ (k, e) = Z- 1 J2 e-K***-^ \ (m\i> a (0)\n)\ 2 

mn 

x {2n) 3 S(k - [P„ - P m ]/H) S(e - [E n - E m \) 

(2.20) 

and 

A-(k,e) = Z- 1 ^]e-' 3 ( B ''-' lAr '')|(m|^(0)|7i)| 2 

mn 

x (27r) 3 <5(k - [P„ - P m ]/h) S(e - [En - E m }) 

(2.21) 

can be expressed in terms of matrix elements of single fermion 
field operators ip a (0) at the origin between the exact eigen- 
states |n) of the many-body system. Here P„, E n , and N n 
are the corresponding eigenvalues of momentum, energy, and 
particle number, respectively. In thermal equilibrium, the par- 
tial spectral functions are related by the detailed balance con- 
dition 



A_(k,e) = e'^-rt A+{k,e) . 



(2.22) 



At zero temperature, therefore, the hole part A- (k, e) of the 
spectral function vanishes for e > (J, and vice versa the particle 
part A + (k, e) vanishes for e < ^. The total spectral weight in 
the hole part 

J deA_(k,e) = n a {k) (2.23) 

at arbitrary temperatures is equal to the fermion occupation 
numbern cr (k) = — G(k,r = —0) for a single spin orientation 
a (in the balanced gas discussed here, both components a = 
±1 have the same occupation, of course). 

Within the BCS description of fermionic superfluids, the 
spectral function consists of two infinitely sharp peaks [40] 



A(k, s) = ui S(e - + vi S(e Ef') 



(2.24) 



which represent the particle and hole part of the spectral func- 
tion. The associated energies 



4 ±} = /i±v> k -M) 2 + A2 



(2.25) 



— ihliJn + £ — [I 



(2.17) 



describe the standard dispersion of the Bogoliubov quasipar- 
ticles. They exhibit a finite gap, whose minimum value A 
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is taken at a finite momentum = ^/2mfi/h (note that 
fi — > £f > in the BCS limit). Within the standard BCS 
theory, these excitations have infinite lifetime at arbitrary mo- 
menta k and there is no broadening or incoherent background. 
Going beyond the exactly solvable BCS Hamiltonian, how- 
ever, gives rise to a finite lifetime of the fermionic excitations 
and thus will broaden the two delta-peaks in Eq. (2.24) even 
at zero temperature. It is our aim in the following, to calculate 
these effects quantitatively for the simple model Hamiltonian 
Eq. (2.1) in the whole range of coupling strengths and temper- 
atures. 

Eq. (2.17) has the form of a Cauchy integral in the theory 
of complex functions. It is therefore convenient to define a 
complex Green function G(k, z) depending on a complex fre- 
quency z, which is analytic in the upper and lower complex 



of the spectral functions A and Af of the initial and final 
states. Here, an unknown prefactor that depends on the in- 
teraction parameters for the coupling to the rf field has been 
set equal to h, which provides a convenient normalization for 
the total weight integrated over all frequencies (see (2.36) be- 
low). This overall constant drops out in normalized spectra 
by dividing out the zeroth moment J duil(q, us) or has - in 
any case - to be adjusted to the measured signal in compar- 
ison with experimental data. Since the wave vector q of the 
rf field is much smaller than those of the atoms, it is an ex- 
cellent approximation to set q = 0. In the absence of a 
probe that selects atoms according to their momenta k, the 
spectrum /(q = 0, w) = I(co) is thus only a function of 
the rf frequency co. In addition, for the standard situation 
with an empty final state /, the partial spectral functions are 
half planes lm(z)<0, respectively. This complex Green func- Af t+ (k,e) = Af (k, e) and (k, e) = 0. Using (2.28), 



tion is related the Matsubara Green function and to the spec- 
tral function by 

Q(k,u n ) = G(k, z = fi/K + iu> n ) , (2.26) 
A(k,e) = ±TT- 1 Im{G(k,z = e/h±iO)} , (2.27) 

respectively. Thus, in a first step we obtain the complex Green 
function G(k, z) as an analytic continuation from the Matsub- 
ara Green function G(k, oj n ). In a second step, we insert the 
complex frequency z = e/h±iO and obtain the spectral func- 
tion A(k, e) from (2.27). The fact that (k, uj n ) uniquely de- 
termines the spectral function has been proven by Baym and 
Mermin [41]. 

In practice the analytic continuation for calculating the 
spectral function A(k, e) is done by using the maximum- 
entropy method [42] which is described in detail in Appendix 
A. We have checked the accuracy of our results a posteriori by 
inserting the calculated spectral functions in Eq. (2.17). The 
given 'initial' data G(k, oj n ) are then found to be reproduced 
with a relative accuracy that is typically in the 10~ 5 range. 



C. Rf response 

In radio-frequency experiments, the external rf field trans- 
fers atoms from one of the two occupied spin states (as initial 
state) into an empty final state. In the following, we assume 
that the final state, which is denoted by an index /, has a neg- 
ligible interaction with the initial one. It can thus be described 
by the free-fermion spectral function 



A f (k > e) = S(e-[E f +e k }) 



(2.28) 



where Ef is the excitation energy of the final state, which 
has a free particle dispersion £k = h 2 k 2 /2m. Within linear 
response, the rate of transitions out of the initial state induced 
by the rf field with frequency ui and wave vector q is given by 
a convolution 

/d 3 k f 
— ^ J de[A f . + (k + q,e + huj)A_(k,e) 

- A ft _(k + q, e + huj)A + {k, e)} 

(2.29) 



the resulting rf spectrum 



d 3 k 
(2^)3 



A_(k, t 



Hut) 



(2.30) 



is an integral over the hole part A_ (k, e) of the single par- 
ticle spectral function in the initial, strongly correlated state. 
For convenience we have taken Ef = 0, which redefines the 
position of zero frequency u = in the rf spectrum. 

Within a BCS description, the spectral function is given by 
Eq. (2.24). Its hole excitation part has a delta-peak at E^ ■* 
whose weight is equal to the occupation number n a (k) = 
This reflects the simple fact that a hole with momentum k can 
only be created if a fermion is present with this momentum. 
The resulting rf spectrum in our normalization is 



-TbcsM 



,3/2 



2 1 /2 7r 2/ l 2 



huj 



A 2 
2hw 



1/2 



A 2 



It exhibits a sharp onset at hw T , 



2{nw) 2 ■ 

(2.31) 
[i. As 



will be shown below, such a sharp onset is not found from our 
numerical results for the spectral function, even in the weak- 
coupling limit v <C — 1. The origin of this discrepancy may be 
traced back to the fact that the dominant contributions to the rf 
spectrum near ui m .i n arise from the spectral function A- (k, e) 
in the limit k — > 0, i.e. far from the Fermi surface at kp. Now, 
deep in the Fermi sea, the true spectral function is not de- 
scribed properly by an extended BCS description, which has 
sharp quasiparticles at arbitrary momenta. In fact, the simple 
form (2.24) of the single fermion spectral function holds only 
if the interaction part of the full Hamiltonian Eq. (2.1) is ap- 
proximated by the exactly soluble reduced BCS Hamiltonian 
[43]. Its interaction term 



H'bcs 



ov ^ 



2V 



(2.32) 



k,k' 



involves only pairs with vanishing total momentum Q = 0. 
This approximation excludes density fluctuations and there- 
fore does not account for the collective Bogoliubov- Anderson 
mode [4]. Moreover, its fermionic quasiparticles are exact 
eigenstates of the reduced BCS Hamiltonian at arbitrary mo- 
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menta. The difference 



c k+Q,cr C ik,-o- c -k'.-cr c k'+Q,cr 



go_ \ - 

a k,k'.Q^O 

(2.33) 

between the full Hamiltonian Eq. (2.1) and that of the reduced 
BCS model therefore gives rise to residual interactions be- 
tween the quasiparticles and their coupling to the collective 
Bogoliubov-Anderson mode. This will be discussed in more 
detail in Appendix B. As will be shown quantitatively in Sec. 
IV, the residual interactions result in an appreciable broad- 
ening 7(k) of the spectral functions. In particular, the hole 
part becomes increasingly broad as k — > (see Fig. 2a for a 
coupling strength v = — 1). A BCS-type rf spectrum (2.31) 
requires that 7(k = 0) <C tui) m i n w A 2 /2ep in weak cou- 
pling. This condition is never fulfilled in practice, because the 
gap A ~ exp (ttv/2) vanishes exponentially in the BCS limit 
v <C — 1, while 7(k = 0) can be shown to be of order (kpa) 2 
as kp\a\ "C 1 due to the decay via intermediate states in- 
volving three quasiparticles (see Eq. (B8) in Appendix B). As 
a result, the onset and peak shift of the rf spectrum in weak 
coupling are dominated by the Hartree contributions and do 
not reflect the appearance of a pairing gap. 

In the limit v ^> +1 of a molecular BEC, the fermions 
form a superfluid of strongly bound dimers. In this regime, 
the gap parameter A becomes negligible compared with the 
magnitude of the chemical potential and the hole excitation 
energy (2.25) approaches ' — > 2/i — e^. Since the ex- 
tended BCS description of the crossover becomes exact again 
in the molecular limit, where it reduces to an ideal Bose gas 
of dimers, one can use Eq. (2.24) for the associated spectral 
function of fermionic excitations, which gives 



A_ (k, e) = vl S(e + e k - 2/i) 



(2.34) 



in the BEC limit. The weight v 2 . = 4nna 3 (l + k 2 a 2 )~ 2 now 
coincides with the square of the bound state wave function in 
momentum space. The resulting rf spectrum 



BEC 



n (Huj + 2/x) 1 / 2 



(2.35) 



is a special case of that derived by Chin and Julienne [44] in 
the molecular limit for bound-free transitions in the absence of 
final state interactions. It has an onset hw m in,BEC = — 2[i — * 
Eb that is determined by the molecular binding energy Eb = 
H 2 /ma 2 , as expected. This energy also sets the scale for the 
half width of the rf spectrum, which is E w =72;, with a 
numerical factor 7 = 1.89. 



D. Rf spectra at high frequencies and contact density 

Our definition of the rf spectrum in (2.30) and the normal- 
ization (2.23) of the hole part of the spectral function imply 
that the total weight integrated over all frequencies 

J duj I(u) = n a = n/2 (2.36) 



is determined by the density n a of atoms from which the 
transfer to the empty final state / occurs. This normaliza- 
tion fixes the overall prefactor and determines the normalized 
form of the rf spectra, in which the zeroth moment is divided 
out. An analysis of the spectra in terms of their nontrivial 
higher moments, however, does not seem to work. Indeed, 
it follows from (2.31) and (2.35) that the rf spectra at high 
frequencies fall off like u; -3 / 2 , both in the BCS and the BEC 
limit. Thus already the first moment of the spectrum diverges. 
The issue of the behavior at high frequencies has been investi- 
gated recently by Schneider et al. [45]. They have shown that 
the exact expression (2.30) quite generally implies an oj~ 3 / 2 
power law decay 



1/2 



-3/2 



(2.37) 



Here, the coefficient C is defined by the behavior n a (k) — > 
C/fc 4 of the momentum distribution at large momenta. It was 
introduced by Tan [23] as a parameter that characterizes quite 
generally fermionic systems with zero range interactions. As 
shown by Braaten and Platter [24], this coefficient is a mea- 
sure of the probability that two fermions with opposite spin 
are close together and is thus called a contact density or sim- 
ply the contact. In the balanced superfluid, it has actually been 
determined from a measurement of the closed channel frac- 
tion by Partridge et al. [46], as analyzed in detail by Werner, 
Tarruell and Castin [47]. In the BEC limit, the well known 
expression for n a (k) in terms of the square of the bound state 
wave function yields Cbec = 47m/a, consistent with the 
explicit form (2.35) of the spectrum in the BEC limit. 

There are two important points in this context, which we 
discuss in the following. First of all, the asymptotic w~ 3 / 2 
power law decay of the exact rf spectrum is valid only in 
the ideal case of zero range interactions and identically van- 
ishing final state effects. Indeed, an explicit calculation of 
the rf spectrum in the molecular limit by Chin and Julienne 
[44] shows that in the presence of a nonzero scattering length 
aj 7^ between the hyperfine state that is not affected by 
the rf pulse and the final state of the rf transition, the spec- 
trum decays like uj~ b / 2 at large frequencies. The short range 
part of the interaction, that is responsible for the slow decay 
of the spectrum, is therefore cancelled out by the interaction 
between the final state and the state that remains after the 
rf transition. This result remains valid quite generally along 
the whole BCS-BEC crossover, as discussed by Zhang and 
Leggett [48, 49]. In particular, this behavior guarantees that 
the rf spectrum has a finite first moment. As shown in Refs. 
[12] and [13], it allows to define an average 'clock shift' 



44 



1 

9f 



(2.38) 



that is again determined by the contact coefficient C = skp 
and the renormalized interaction constants g = 4irh 2 a/m and 
gf = 4irh 2 af/m. In particular, there is a perfect 'atomic 
peak' I{u>) ~ S(u>) and no clock shift at all if a = a,f. The 
existence of higher moments of the rf spectrum relies on ac- 
counting for the nonzero range ro 7^ of the interaction. 
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Since this is expected to affect the spectrum only at frequen- 
cies of order h/mr 2 , this regime, however, will hardly be ac- 
cessible experimentally. 

As a second point, we consider the behavior of the contact 
coefficient C in the weak coupling limit. Standard BCS theory 
for the momentum distribution at large k predicts that the cor- 
responding dimensionless factor s = C/k F is exponentially 
small sbcs = (A/2ef) 2 - This is in agreement with the high 
frequency asymptotics of the ideal BCS spectrum (2.3 1) with- 
out final state interactions according to the result in (2.37). 

It turns out, however, that the exponentially small value of 
the contact coefficient in the BCS limit is an artefact of work- 
ing with a reduced BCS Hamiltonian, which only takes into 
account the pairing part (2.32) of the interaction. By contrast, 
the full Hamiltonian gives an additional contribution that is as- 
sociated with non-condensed close pairs. For weak coupling, 
this is much larger than that of the condensed pairs described 
by the BCS theory. More precisely, it turns out that the coef- 
ficient 

s = [A 2 - ru (r = 0, r = -0)]/(4eJ.) (2.39) 

in front of the n a (k) = s(fci?/fc) 4 behavior of the momen- 
tum distribution at large k contains a nontrivial contribution 
associated with the upper diagonal element Tn of the ver- 
tex function defined in (2.12) in the limit of vanishing spa- 
tial and temporal separation. In the molecular limit, the con- 
tribution from this term is negligible. The asymptotic re- 
sult Abec = 2e^-\/4u/37r for the gap parameter then gives 
rise to a linearly increasing dimensionless contact parameter 
sb ec = 4w/37r, consistent with the naive result discussed 
above. On the contrary, in the weak coupling limit, the con- 
tribution from non-condensed close pairs is dominant com- 
pared with the exponentially small BCS contribution from 
condensed pairs. In the limit v <C — 1 the leading behavior 
is given by 

-Tii (0,-0)= . (2.40) 

The resulting dimensionless contact coefficient s wc = 
(2/3irv) 2 in weak coupling is therefore much larger than the 
exponentially small BCS contribution. It is remarkable that 
the leading term in the weak coupling contact density of the 
superfiuid with a < is identical to the one that is obtained 
in a repulsive dilute normal Fermi liquid with a > 0, that has 
first been calculated by Belyakov [50]. This shows that the 
dominant contribution to the contact density is independent of 
the sign of the interaction, consistent with the 'adiabatic theo- 
rem' 

am = <"'> 

that relates the derivative of the energy per volume u with re- 
spect to the inverse scattering length to the contact coefficient 
C [24, 51]. In fact, the simple mean-field interaction energy 
linear in a, which is the leading correction to the ground state 
energy of the ideal Fermi gas, shows that C(a) ~ a 2 is inde- 
pendent of the sign of interactions to lowest order. The BCS 




V = l/k„a 

F 

FIG. 1: (Color online) The dimensionless contact coefficient s is 
shown as a function of the dimensionless coupling strength v = 
1/k.Fa. The red solid line represents our numerical result obtained 
from Eq. (2.39). The left and right black dashed lines represent the 
asymptotic formulas s wc and sbec given in the text, respectively. 



pairing effects, that appear in the case of a negative scatter- 
ing length, only give a subdominant, exponentially small re- 
duction of the energy that is reflected in a corresponding tiny 
enhancement of the contact density. The full dependence of 
s(v) along the BCS to BEC crossover for the balanced gas 
at zero temperature is shown in Fig. 1. The particular value 
s(0) = 0.098 at unitarity has in fact been determined before in 
the context of the average clock shift (2.38) [12] and is close 
to our present value s(0) = 0.102 that follows from (2.39). 
Since the contact density is a short range correlation property, 
it is not very sensitive to temperature, see [52]. 

An important consequence of the failure of naive BCS the- 
ory to account for the correct value of the contact density C 
in weak coupling is the fact that the weight of the rf spectrum 
at high frequencies that is determined by Eq. (2.37) in the ab- 
sence of final state interactions or by the average clock shift 
(2.38) is strongly underestimated by using the idealized form 
(2.24) of the spectral functions that follow from a naive BCS 
theory. It is an interesting open problem to determine analyt- 
ically the explicit form of the spectral functions in the weak 
coupling limit which is consistent with the correct high fre- 
quency asymptotics (2.37) with the proper value for the con- 
tact density. 



III. EFFECTIVE FIELD THEORY AND PAIR SIZE 

For a molecular BEC that consists of tightly bound dimers, 
the notion of a pair size is well defined. In the relevant case of 
a zero range two-particle interaction with positive scattering 
length a, it is determined by the rms extension £ m = a/v2 
of the two-body bound state. Since kpa <C 1 in the BEC 
limit, the size of the pairs in this regime is much smaller 
than the average interparticle spacing (37T 2 ) 1 / 3 /kp w 3.1/fcp 
of the fermionic gas in the absence of an attractive inter- 
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action. Motivated by the fact that rf spectroscopy in this 
limit effectively reduces to a two-body molecular spectrum, 
a spectroscopic pair size £ w has been determined from the 
half width of the measured rf spectrum I(u>) by the relation 
£w = 7 x h 2 /2mE w [18]. By its definition, this pair size co- 
incides with the molecular size £ m = aj\[2 in the appropriate 
limit. Extending this definition to arbitrary coupling strengths, 
one obtains a ground state pair size £ w = 0.86Hvf/ &o in the 
opposite BCS limit, which correctly describes the exponen- 
tially large size of Cooper pairs characteristic for weak cou- 
pling BCS superconductors [53]. It is thus plausible to use 
this spectroscopic definition of the pair size for the complete 
range of couplings, in particular also near the unitarity limit 
[18]. The measured rf spectrum at the lowest temperature, 
that has been reached experimentally, is then found to give 
an effective pair size £ w ~ 2.6/kp at unitarity [18]. This is 
somewhat smaller than the average interparticle spacing, in- 
dicating that the unitary gas has pairs that no longer overlap. 
The numerical value for the pair size is in fact close to that 
obtained from calculating the width of the pair wave function 
in a variational Ansatz for the ground state of the BCS-BEC 
crossover [54]. 

An obvious question in this context is, whether there are in- 
dependent measures of the pair size, which do not rely on the 
spectroscopic definition that is motivated by the extrapolation 
from the molecular limit or on the related variational approxi- 
mation for the ground state in terms of a product of two-body 
wave functions. In the following, we will show that a many- 
body definition of the pair size can be obtained from the q- 
dependent superfluid response function, following the basic 
concept of a nonlocal penetration depth in superconductors 
[53]. This response can be calculated from an effective field 
theory of the superfluid state, including the next-to-leading 
order corrections to the standard quantum hydrodynamic La- 
grangian. Remarkably, the value of the pair size at unitarity 
that follows from the g-dependent superfluid response is close 
to that inferred spectroscopically from the half width of the rf 
spectrum. 

The basic idea, that allows to define a characteristic length 
£ p of a fermionic superfluid without reference to an ap- 
proximate BCS or molecular description of the many-body 
ground state is related to the well known calculation of the 
q-dependent penetration depth A(g) in charged superconduc- 
tors. The latter is defined by the nonlocal generalization 
j(q) = — A(q)/4irA 2 (q) of the London equation, relating the 
current density induced by a transverse vector potential in lin- 
ear response [53]. The square of the inverse effective pene- 
tration depth determines the superfluid density n s (q), which 
obeys n s (q = 0) = n in any Galilei invariant superfluid. For 
finite momentum q the superfluid response is reduced by a 
correction that has to vanish like q 2 in an isotropic system. 
The correction defines a characteristic length £ p according to 



n s (q) 



" 2 >-9 

30^ 



(3.1) 



in the limit of small wave vectors q. Here, the prefactor 
in the q 2 -correction has been chosen in such a way, that 



the characteristic length £ p coincides with the Pippard length 
£p = Hvp/nAo in the weak coupling limit . 

In order to determine the value of £ p at unitarity, one needs 
the leading order corrections in an expansion in small gradi- 
ents to the universal quantum hydrodynamic Lagrangian den- 
sity 



Cn = 



2m 



V-(v?) s 



(3.2) 



of a translation invariant, neutral superfluid with (Bogoliubov- 
Anderson) sound velocity c s . For the unitary Fermi gas, where 
c 2 = 2/i/3m exactly, these corrections have been discussed 
in detail by Son and Wingate [25]. Restricting ourselves to 
the harmonic description of the Goldstone mode described by 
(3.2) to leading order, the next-to-leading terms in the effec- 
tive field theory are of the form [25] 



£ 



C2 



m 



(3.3) 



The associated dimensionless coefficients ci,2 can only be de- 
termined from a microscopic theory. Their physical meaning 
becomes evident from the fact that c-i determines the reduc- 
tion of the superfluid response for finite wave vectors q as de- 
scribed in (3. 1). In terms of the pair size £ p defined there, one 
finds 



3b^ = 9C2 " 



hn 



(3.4) 



which also makes clear that C2 has to be positive. In contrast 
to C2, the coefficient c\ has no direct physical interpretation. 
From the plane-wave solution of the linear equations of mo- 
tion for the phase fluctuations that follow from the total La- 
grangian density Co + C it is easy to see, however, that this 
coefficient appears in the next-to-leading corrections in the 
dispersion uj{q) = c s q(l — aq 2 jk 2 F + . . .) of the Bogoliubov- 
Anderson mode with a dimensionless coefficient [25] 



;C 2 + Ci 



(3.5) 



Here £(0) ~ 0.36 is the Bertsch parameter, that relates the 
sound and bare Fermi velocities by c 2 = £(0)v F /3. Simi- 
lar to the Bertsch parameter, which appears in the leading or- 
der Lagrangian (3.2), the coefficients ci t 2 can be calculated in 
an expansion around the upper critical dimension four of the 
unitary Fermi gas, as suggested originally by Nussinov and 
Nussinov [55] and started by Nishida and Son [56]. A one 
loop calculation of the coefficients ci$ has recently been per- 
formed by Rupak and Schafer [57]. The resulting value of c\ 
at e = 4 — d = 1 turns out to be c\ f» —0.02. Unfortunately, 
for C2, the one-loop calculation is not sufficient, because a fi- 
nite value of C2 only appears at order e 2 [57]. This is easy to 
understand from the connection (3.4) between C2 and the pair 
size, which is expected to vanish linearly in e = 4 — d. Indeed 
in four dimensions, a two-particle bound state in a zero range 
potential only appears at infinitely strong attraction [55]. The 
unitary gas in d — 4 therefore has a vanishing dimer size and 
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is effectively an ideal Bose gas, similar to the situation in the 
BEC limit in d = 3 [58]. In order to fix the value of C2 for 
the unitary gas in three-dimensions, we use the connection 
(3.5) between the next-to-leading order coefficients of the ef- 
fective field theory and the ^-corrections to the dispersion 
u>(q) of the Bogoliubov- Anderson mode. This dispersion has 
been calculated within a Gaussian fluctuation approximation 
for arbitrary coupling strengths v [39] and exhibits a negative 
curvature with a « 0.06 right at unitarity [59]. Combined 
with the value of c\ from the e-expansion, this leads to the es- 
timate C2 ~ 0.02 for the unitary gas in three dimensions. As 
a result, the pair size that follows from (3.4) turns out to be 
£ p « 2.Q2/kp. It is remarkable that this value essentially co- 
incides with that inferred from the spectroscopic definition in 
Ref. [18] or the width of the pair wave function in Ref. [54]. 
It should be noted, though, that apart from the uncertainties in 
the precise values of ci^, there is a certain amount of arbitrari- 
ness in defining a 'pair size' , both from the rf spectrum or from 
the g-dependent superfluid response via (3.4). This is related 
to the precise choice of the prefactor both in the spectroscopic 
definition and in (3.4), where the Pippard length has been used 
as a reference scale. The naive conclusion that the unitary gas 
has non-overlapping pairs in its ground state should therefore 
be viewed with a great deal of caution. 

The fact that the coefficients ci 2 of the next-to-leading or- 
der Lagrangian (3.3) have a comparable magnitude at unitar- 
ity has a further interesting consequence. Indeed, these coef- 
ficients determine the magnitude of the Weizsacker inhomo- 
geneity correction [60, 61] 



; H 2 (Vn(x)Y 
2m n(x) 



(3.6) 



to the ground state energy density of the unitary Fermi gas 
[57]. The associated dimensionless coefficient b is related to 
the q 2 -corrections of the density response [25] 



X (q) = X(0) 1 - b 



hq 



For the unitary Fermi gas, the coefficient b is again determined 
by the next-to-leading order Lagrangian (3.4) via [57] 



Cl 



(3.8) 



(3.7) 



Using our estimates for 01,2, this leads to b f=a 0.22, a value 
that is almost one order of magnitude larger than the result 
&(°) = 1/36 = 0.028 which is obtained for an ideal Fermi 
gas. The unitary gas is therefore remarkable in the sense that 
its kinetic energy density Eki n {x) ~ £(0) n 5 ^ 3 (x) is reduced 
by the Bertsch parameter £(0) « 0.36 compared with the non- 
interacting gas, yet the coefficient b of the Weizsacker inho- 
mogeneity correction is strongly enhanced [57]. 



IV. NUMERICAL RESULTS 

In the following, we present numerical results for the spec- 
tral function A(k, e) and the rf spectrum I{u>) which can be 
compared with experimental data. The numerical calcula- 
tions are performed in two steps. First, the Matsubara Green 
function (k, u n ) is calculated by solving the self-consistent 
equations (2.10)-(2.15). In a second step the spectral func- 
tion A(k, e) is calculated from (2.17) by analytical continu- 
ation as described in Sec. II B. For this purpose we employ 
a maximum-entropy method that is described in Appendix A. 
Eventually, the rf spectrum I(uj) is calculated by evaluating 
the momentum integral (2.30) numerically. 



A. Spectral functions 

Our numerical results for the spectral functions A(k, e) 
in the experimentally relevant range of interaction strengths 
v = —1,0 and + 1 are shown in Figure 2. The associated 
temperature is T = 0.01 Tp, i.e. deep in the superfluid regime 
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FIG. 3: (Color) Density plots of the spectral function A(k, e) at unitarity (v — l/(kFd) = 0) for different temperatures. From top left to 
bottom right: T/T F = 0.01, 0.06, 0.14, 0.160(T C ), 0.18, 0.30. The white horizontal lines mark the chemical potential fj,. At temperatures 
smaller than the superfiuid transition temperature T c two quasiparticle structures with a BCS-like dispersion can be seen. The width of the 
spectral peaks is of the same order as the quasiparticle energy. With increasing temperature the two branches gradually merge into a single 
quasiparticle structure with a quadratic dispersion above T c . Note however, that the quadratic dispersion is shifted to negative frequencies 
compared to the bare Fermion dispersion relation. This Hartree shift is of the order of U — — 0.46 £f and is essentially responsible for the 
shifted rf spectra in the normal phase in Fig. 6. 



in all three cases. Evidently, both at v = —1 and at unitar- 
ity v = 0, a BCS-like quasiparticle structure appears, with an 
excitation gap whose minimum is at a finite value of the mo- 
mentum. On the BEC side, at v = +1, the backbending in 
the dispersion curve has apparently disappeared. This is con- 
sistent with the expected existence of a critical value v s > 0, 
beyond which the fermionic excitations have their minimum 
at k = 0. From our numerical data on the momentum de- 
pendence of the fermionic excitation spectrum, the associated 
critical coupling constant is v s = 0.8. This is about a factor 
of two larger than the mean-field prediction, which is deter- 
mined by the zero crossing of the chemical potential. The fact 
that v s occurs in the regime where the chemical potential is 
already negative has been noted before in an e = 4 — d expan- 
sion by Nishida and Son [62]. Extrapolating their one loop 
result to e = 1 gives fj, s ~ — 0.5 ef at the critical coupling 
v s , in rather good agreement with our result [i a = — 0.5Aef- 
In population imbalanced gases the change in the curvature of 
the fermionic excitation spectrum at v s determines the criti- 
cal coupling of the splitting point S, at which the continuous 
transition from a balanced to a an imbalanced superfiuid on 
the BEC side splits into two first order transitions [63]. 









particle 


hole 


V 




A/e F 


m*/m U/ef 


77i* /m U/ef 


-1 


0.73 


0.14 


1.05 -0.26 


1.12 -0.17 





0.36 


0.46 


1.00 -0.50 


1.19 -0.35 


+1 


-0.93 


1.10 


1.02 -0.42 


1.28 -0.37 



TABLE I: Effective mass m* and Hartree shift U of the quasiparticle 
dispersion relations at T — 0.01 Tf, obtained by fitting Eq. (4.1) to 
the peak maxima of the spectral functions in Fig. 2. 

Empirically, the form of the quasiparticle dispersion rela- 
tions may be extracted from the peak position of the spectral 
function. It turns out that these peaks fit reasonably well to a 
modified dispersion 

^ ±) =M±^(^k + tf-/,) a + A* (4.D 

of Bogoliubov quasiparticles, in which the effective mass m* 
and an additional Hartree shift U are used as fit parameters. 
The associated values for m* and U that follow from the spec- 
tral functions shown in Fig. 2 are summarized in Tab. I. It is 
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interesting to note that both the masses and the Hartree shifts 
are different for particle and hole excitations. The particle- 
hole symmetry of the standard BCS description of the quasi- 
particle dispersion is therefore broken at these large coupling 
strengths. A second feature of interest is that the hole disper- 
sion relation starts to deviate from the BCS form (4.1) 
only for momenta k > 1.5 kF, when the spectral weight of 
the hole peak is smaller than 0.5%. 

Using QMC methods, the particle dispersion relation at uni- 
tarity and T = has been calculated previously by Carlson et 
al. [64]. Our values for the Hartree shift U and the effec- 
tive mass m* of the particle dispersion relation agree reason- 
ably well with the QMC values. Experimentally, the Hartree 
shift U was extracted recently from rf measurements by Schi- 
rotzek et al. [65]. In this work, the measured peak positions 
of the rf spectra were fitted with the peak position obtained 
from the BCS formula (2.31), including an additional Hartree 
shift: = 4( v /( A i - U) 2 + 15A 2 /16-/^ + [//4)/3. If we 
apply this method to our calculated rf spectra, we obtain dif- 
ferent values for U than those listed in Tab. I. In particular, this 
method gives U = -0.28, -0.52, -0.22 for v = -1, 0, +1 
at T = 0.01 Tp and doesn't take the effective mass into ac- 
count (we note, that the rf spectrum is only sensitive to the 
hole excitation part of the spectral function). This discrep- 
ancy is probably due to the fact, that the assumption of having 
sharp quasiparticles is not reliable in this regime. 

It is evident from the quantitative form of the spectral func- 
tions, that the parametrization of the fermionic excitations by 
a modified dispersion (4.1) is not an adequate description of 
the excitation spectrum because of the rather strong broaden- 
ing of the quasiparticle peaks, even at very low temperatures. 
The physical origin of this broadening is the residual interac- 
tion between the quasiparticles that follows from the Hamil- 
tonian in Eq. (2.33). As shown in Appendix B, this interac- 
tion leads to a finite width of the spectral functions even at 
T = 0, except near the minimum of the dispersion curve for 
the particle excitations and close to the maximum of the dis- 
persion for the hole excitations. Here, the quasiparticle life- 
time broadening has to vanish because there are no available 
final states into which it may decay. Focusing on the inter- 
action of the quasiparticles with the collective Bogoliubov- 
Anderson mode, which is the dominant mechanism for decay 
near the minimum (maximum) of the particle (hole) disper- 
sion curve, it is straightforward to see that there is actually a 
finite interval in momentum space where the spectral width 
vanishes identically. This width is determined by the kine- 
matic constraint that that the quasiparticle decay by emission 
of phonons is possible only if the group velocity dE^/Ok of 
the fermionic excitations becomes larger than the sound ve- 
locity c s . In fact a similar situation appears for a hole in 
a Neel ordered antiferromagnet, whose spectral function is 
sharp as long as its group velocity is below the spin wave ve- 
locity of antiferromagnetic magnons [66]. The fact that our 
numerically calculated spectral functions A(k, e) exhibit a fi- 
nite broadening at the lowest temperature T = Q.QlTp even 
near the dispersion minimum is probably related both to the 
numerical procedure of evaluating A(k, e) using the Maxent 
technique which can never give rise to perfectly sharp peaks, 




FIG. 4: (Color online) The spectral function A(k, e) as a function 
of e for selected fixed values k at unitarity v = l/(/cpa) = and 
at criticality T/Tf = 0.160(T C ). The selected values of the wave 
number k are represented by the colors of the lines corresponding to 
the peaks from left to right: k/k F = 0.00 (black), 0.52 (red), 0.77 
(orange), 1.00 (green), 1.26 (cyan), 1.51 (blue), 2.02 (magenta). 
The different methods for calculating the spectral function are dis- 
tinguished by the line styles: maximum-entropy method (solid lines) 
and Pade approximation (dashed lines). 



but also to the self-consistent structure of our Luttinger-Ward 
formulation. Indeed, in a diagrammatic language, the latter 
implies summation of diagrams with identical intermediate 
states for Fermions, which - in an exact theory - are excluded 
by the Pauli principle. Unfortunately, to our knowledge, there 
exist no analytical results on the broadening of the Bogoliubov 
quasiparticles beyond the perturbative treatment outlined in 
Appendix B. Experimentally, this question may in principle 
be resolved by studying momentum resolved rf spectra that 
have recently been obtained by Stewart et al. [19]. Unfortu- 
nately, at present, experimental data on spectral functions are 
available only near the critical temperature of the superfluid 
transition, where the finite lifetime arises due to the scattering 
with thermally excited quasiparticles. 

To discuss the situation at finite temperature, we plot the 
spectral function A(k, e) at unitarity for different tempera- 
tures above and below T c in Figs. 3 and 4. It is interesting 
to observe how the two BCS-like quasiparticle peaks evolve 
with increasing temperature and finally merge into a single ex- 
citation structure with a quadratic dispersion at temperatures 
around T c . Note however, that the spectral peak in the nor- 
mal phase is shifted to negative energies compared to the free 
Fermion dispersion relation = ft 2 k 2 /2m. This Hartree 
shift is responsible for the observation of shifts in experimen- 
tally measured rf spectra above T c , that will be discussed in 
detail below. The observation of such a shift in the rf spectra 
in the normal state is therefore not necessarily a signature of 
pseudogap effects. 

Finite temperature QMC calculations of the spectral func- 
tion at unitarity by Bulgac et al. [67] indicate the presence of 
a gapped particle excitation spectrum of the form (4.1) also 
above the critical temperature, which is not found in our ap- 
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proach. More generally, it is evident from the spectral func- 
tions of the unitary gas above T c which are shown in Fig. 3, 
that a simple pseudogap Ansatz for the spectral function [69] 
is not consistent with our results. As can be seen from the 
lower three graphs in Fig. 3, our approach leads to a single, 
broad, ungapped excitation peak with a quadratic dispersion at 
temperatures T >T C instead of two excitation branches with a 
gapped, BCS-like dispersion, as expected from the pseudogap 
approach. In particular we do not observe a strong suppres- 
sion of spectral weight near the chemical potential. 

Apart from the dominant peaks discussed above our spec- 
tral functions show some additional structure that have much 
smaller weight, however. Specifically, at unitarity and tem- 
peratures above T c a small second peak is visible for k < kp 
in Fig. 3. At T = 0.3 Tp this residual peak contains ~ 17% 
of the spectral weight. The situation is similar on the BEC 
side of the Feshbach-resonance at v = 1, where above T c a 
second peak at negative energies is present for k < kp, with 
a spectral weight of ~ 22%. 

Recent experiments by Stewart et al. [19] have succeeded 
to perform rf spectroscopy in a momentum resolved manner, 
from which one directly obtains the hole spectral function 
-A-(k, e) as a function of both, momentum and energy. A 
quantitative comparison with our calculated spectral functions 
is difficult however, since the measured spectral functions in- 
volve an average over the inhomogeneous density profile of 
the trapped atoms. Nevertheless, as shown in Fig. 5, the qual- 
itative structure of our hole spectral function of the uniform 
system near the critical temperature is similar to that observed 
experimentally. To separate the intrinsic from an inhomoge- 
neous broadening in a trap requires to combine momentum 
and local resolution, which is currently investigated in the 
group at JILA. Experiments of this kind would allow to dis- 
tinguish between different models for the spectral functions, 
in particular for the 'pseudogap' phase immediately above T c . 
The existence of preformed pairs in this regime is often de- 
scribed by sharp spectral functions of the form (2.24) with a 
nonvanishing gap parameter A pg . As shown recently by Chen 
et al. [70] and Dao et al. [71], this assumption is also con- 
sistent with present experimental data, due to the inhomoge- 
neous averaging associated with the position dependent gap 
parameter in a trap. 



B. Rf response 

In Fig. 6 we show the calculated rf spectra, together with 
the locally resolved experimental data of the MIT group [65]. 
The measured rf data shown in Fig. 6 have been corrected for 
the small mean-field final state interaction energy, which al- 
lows for a direct comparison with our calculated spectra. For 
a detailed comparison we must take the finite experimental 
resolution into account, however. The MIT group uses an ap- 
proximately rectangular rf pulse with a length of T = 200/is 
in order to transfer atoms to the empty hyperfine state. Thus, 
the Fourier spectrum of the radio-frequency source has a fi- 
nite width and the calculated spectra obtained using Eq. (2.30) 
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FIG. 5: (Color) Density plot of the hole part of the spectral function 
A_(k,e) at unitarity v = l/{k F a) = mdT/T F = 0.150 slightly 
below T c (left) and in the BEC-regime at v = l/(kpa) = +1 
and T/Tf — 0.207 at the superfluid transition temperature (right). 
The white horizontal line marks the chemical potential fj,. The color 
scheme is the same as in Figs. 2 and 3. 

need to be convolved with sinc 2 (wT/2), i.e. 

7 e xp(w) = J dcu' I(u - w') sinc 2 (w'T/2) . (4.2) 

The finite experimental resolution thus leads to a slight broad- 
ening and a small shift to higher frequencies of the calculated 
rf spectra. At unitarity and T = 0.017V, the broadening is 
- 0.07 E F and the shift - 0.01 e F . 

As can be seen in Fig. 6, the rf spectra in the homogeneous 
system show a single peak that is shifted compared to the bare 
transition frequency, which is set to uj = for convenience. 
Apart from slightly overestimating the width of the spectral 
lines, our numerically obtained spectra agree very well with 
the experimental data. Note that no fitting parameters have 
been used, apart from adjusting the absolute height. 

In the first rf measurements by Chin et al. [10] a secondary 
peak at the bare transition frequency has been observed and 
attributed to the presence of unpaired atoms. In these exper- 
iments however, the measurement of the spectra involved an 
average over the inhomogeneous density profile of the trapped 
atoms. Since more recent locally resolved rf measurements 
[17] didn't show signs of an atomic peak, it is likely that these 
peaks either originate from the low density regions at the edge 
of the atomic cloud, or are an effect of the strong final state 
interactions. 

It is important to notice that there are essentially two con- 
tributions to the rf peak shift. The first one is due to pair- 
ing correlations, which are particularly important in the su- 
perfluid phase and give the dominant contribution to the peak 
shift on the BEC side of the crossover, where the Fermions 
are paired in two-body bound states. The second contribution 
comes from Hartree-type correlations. The Hartree contribu- 
tion dominates the peak shift in the normal phase, where pair- 
ing correlations are small. Furthermore, it also dominates in 
the BCS regime, since the gap is exponentially small whereas 
the Hartree contribution scales linearly with kpa. Together 
with the discussion of the residual interaction between the 
quasiparticles in section II C this implies, that the BCS for- 
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FIG. 6: (Color online) Comparison of the calculated rf spectra at unitarity with the experimental data of the MIT group [65] at different 
temperatures T. Our numerical result is shown by the red solid line. The experimental data are represented by open diamonds connected by 
straight thin dashed lines. Apart from adjusting the peak heights, no fitting parameters have been used. 



mula (2.31) for the rf spectrum in weak coupling is completely 
misleading. 



V. DISCUSSION 

The results presented above on the spectral functions and 
the associated rf spectra of ultracold fermionic gases near 
a Feshbach resonance have two major aspects. First of all, 
they provide a quantitative description of recent experiments 
on the excitation spectra of the near unitary gas. The the- 
ory covers the complete range of temperatures from the su- 
perfiuid near zero temperature to the anomalous normal state 
above T c which is characterized by strong pairing fluctuations. 
From a theory point of view, our results are of relevance as 
a simple example, where the standard quasiparticle descrip- 
tion of BCS theory is strongly modified. As a result of inter- 
actions between quasiparticles and their coupling to the col- 
lective Bogoliubov-Anderson phonon, the spectral functions 
aquire a finite broadening even at zero temperature, except 
for a small range of momenta around the dispersion minimum 
(maximum) of the particle (hole) excitations. Effectively, the 
fermionic excitations along the BCS-BEC crossover are not a 
Fermi gas as in BCS theory, but are described by a Fermi liq- 
uid picture. The spectral functions therefore have vanishing 
width at zero temperature at a sharply defined surface in mo- 
mentum space, where the excitation energy has its minimum 
[68]. Within a perturbative calculation around the BCS limit, 
this has been shown explicitely in Appendix B. More gener- 
ally, it is expected to be valid for arbitrary coupling because 



the phase space for quasiparticle decay vanishes near the dis- 
persion minimum. We are not aware however, of a general 
proof of this statement for the simple model (2.1) of attrac- 
tively interacting fermions studied here. 

As discussed in Sec. II C, the finite lifetime of the fermionic 
excitations at small momenta is particularly important for the 
onset of the rf spectra shown in Fig. 6. In the BCS picture, 
where the spectral function has vanishing width at arbitrary 
momentum k, the rf spectra would exhibit a sharp onset. As 
argued above, however, the fermionic excitations near k = 
have a finite width even at T = up to the critical coupling 
v s , because they are far away from the maximum of the hole 
dispersion. As a result, the rf spectra have no sharp onset, in 
agreement with the experimental observation. A further im- 
portant aspect of our results is that the naive description of 
the BCS-BEC crossover problem by an extended BCS Ansatz 
[4], that appears to work qualitatively at least for the ground 
state is completely inadequate as far as dynamical correlations 
are concerned. In particular, the simple form of the spectral 
function in Eq. (2.24) that follows from a naive BCS theory is 
never valid because the pure pairing Hamiltonian on which it 
is based misses both the broadening e.g. due to collective ex- 
citations and the large contribution to the contact coefficient 
due to non-condensed close pairs. A rather surprising conclu- 
sion of our work is that the next-to-leading terms in the effec- 
tive field theory for the Bogoliubov-Anderson mode allow to 
give a many -body definition of the pair size which agrees quite 
well with the result found experimentally from the half-width 
of the rf spectrum. 

There are of course a number of open problems which 
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should be addressed in future work. In particular, it would 
be interesting to understand to which extent the normal phase 
above T c can be understood in terms of a pseudogap model, 
which has been applied with reasonable success to understand 
ARPES experiments in the context of high T c cuprates [69]. 
As far as pseudogap-effects are concerned, our spectral func- 
tions close to T c show a different behavior than previous non- 
selfconsistent calculations [72, 73, 74] and more recent QMC 
calculations [67], which exhibit a pseudogap. Quite generally 
it is known, that non-selfconsistent calculations favor pseu- 
dogap behavior, whereas selfconsistent calculations suppress 
pseudogap effects (c.f. [75] and references therein). We em- 
phasize however, that in the context of ultracold gases, the 
available momentum and energy resolution in experiments at 
present is not good enough to map out the spectral function 
in sufficient detail. Apart from the rather good agreement 
between the observed rf spectra and our results for the un- 
derlying spectral function, the confidence that our selfcon- 
sistent Luttinger-Ward approach to the BCS-BEC crossover 
gives quantitatively reliable results is supported by the very 
precise description it provides for thermodynamic properties 
(see the discussion at the end of Section II A), much better 
than non-selfconsistent approaches. It is an open problem to 
determine spectral functions e.g. from QMC data, at the level 
of accuracy that has now been achieved for equilibrium prop- 
erties. 



linear equations which can be solved by standard numerical 
methods. 

Unfortunately, the linear equations are nearly singular even 
if the number of discrete energies e (number of unknown vari- 
ables) is smaller than l max (number of equations). Many 
eigenvalues of the linear matrix are very close to zero, so 
that small numerical errors in the Matsubara Green function 
Q (k, oj n ) are enhanced exponentially. As a result, the spectral 
function A(k, e) can not be calculated by this simple method. 

In order to improve and stabilize the method we need some 
prior information. We assume that a smooth background-like 
prior spectral function Aq{e) is given. We define the entropy 

S(k) = Jde [A(k, e) - A (e) - A(k, e) ln[A(k, e)/A (e)}} 

(Al) 

where the integral is evaluated numerically by the trapezoid 
formula for the above defined discrete energies e. The spec- 
tral function A(k, e) is obtained by maximizing the entropy 
(Al) for given wave vectors k where Eq. (2.17) is used as a 
constraint. This procedure is known as the maximum-entropy 
method and can be derived by Bayes inference [76]. It has 
been applied successfully for calculating the spectral func- 
tion j4(k, e) from the Matsubara Green function G(k, u> n ) in 
Monte Carlo simulations [42]. 

In order to implement the constraints we define the chi 
square 
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APPENDIX A: MAXIMUM-ENTROPY METHOD 

In order to solve Eq. (2.17) explicitly for A(k,e) the in- 
tegral is discretized. We chose equally spaced energies e in 
the inner interval — 10 £p < e < +10 £p and logarithmically 
spaced energies in the outer regions — 10 6 £p < e < — 10 £p 
and +10 Bp < £ < +10 6 Ep. We evaluate the integral by us- 
ing the trapezoid formula. On the right-hand side of Eq. (2. 1 7) 
the Matsubara Green function G(k, ui n ) is given for selected 
Matsubara frequencies a4 on a logarithmic scale, which are 
ordered according to < u>n < < ■ ■ ■ < wi ~ 
10 6 £f- In this way, Eq. (2.17) is transformed into a set of 
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A(k,e) 
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(A3) 

is a dimensionless difference and a is a dimensionless stan- 
dard deviation. Minimizing [x(k)] 2 we recover the constraint 
equations (2.17). 

By solving the self-consistent equations of Sec. II we cal- 
culate the Matsubara Green function Q (k, u> n ) with a rela- 
tive accuracy of about 10~ 5 . For this reason, we expect 
|cZ(k, o;„)| ~ 10~ 5 and chose the fixed value a = 10~ 5 for 
the standard deviation. As a result we observe x(k) ~ 1 in 
our numerical calculations where variations occur by a factor 
of 10 for different wave vectors k. 

Using Bayes inference [76] it can be shown that 



1 



Q(k)=aS(k)--[x(k)] 



(A4) 



is the functional which must be maximized by variation of the 
spectral function A(k, e) for every fixed wave vector k. The 
related necessary condition is 



6Q{k) 
5A(k, e) 



= 



(A5) 



which implies the equations to be solved numerically for 
A(k, e). In Eq. (A4) a is a Lagrange parameter which bal- 
ances the weight between the entropy (Al) and the constraints 
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(A2). For a = we recover the constraint equations (2.17). 
On the other hand, in the limit a — > oo we obtain the prior 
spectral function A(k, e) = Aq(e). Thus, a is a parameter 
which must be adjusted to an intermediate value in order to 
obtain an optimum result for the spectral function A(k, e). 
For low values a the constraints are overweighted. A more 
accurate result is obtained for A(k, e), however, instabilities 
may occur. On the other hand, for higher values a the en- 
tropy is overweighted. A more stable result is obtained which 
however, may be less accurate. 

We have defined the entropy (Al) and the chi square (A2) 
as dimensionless quantities which are of order unity. For this 
reason, we expect that a must be of order unity, too. Actually, 
we find that a = 1 is an optimum choice in most areas of the 
phase diagram except for low temperatures. For this reason, 
we use a = 1 in most cases. However, for low temperatures 
T < 0.5 T c the spectral function A(k, e) is very close to zero 
in the gap region. Since the numerical algorithm considers 
the logarithm ln[A(k, e)] an instability occurs. Hence, in this 
latter case for low temperatures we choose a = 100 in the 
crossover and BEC regime, and a = 1000 in the BCS regime. 

For the success of the method an appropriate choice for the 
prior spectrum Aq (e) is very important. First of all the prior 
spectrum Aq(e) should be a smooth function of the energy 
e which models a broad background spectrum. The special 
form of the entropy (Al) does not require Aq(e) to be nor- 
malized. The constraint equations (2.17) will determine the 
spectral function A(k, e) for small and intermediate energies 
in the inner interval — 10 £p % £ ^ +10 £f- However, the 
constraints will provide less information in the tail regions 
e <C — 10 £f and e ^S> +10 ep. For this reason our method 
is considerably improved if the prior spectrum Ao(e) already 
shows the correct wings for e — ► ±oo. 

Investigating the Matsubara Green function for large Mat- 
subara frequencies oj n — ► ±oo we obtain the asymptotic for- 
mula 

C/(k, Lu n ) ?s (— ihwn)^ 1 + a'k(-'ifi^n) 2 + irb(—ihui n )~ 5 / 2 

(A6) 

where a k = -(e k - n) and b = (2/3tt 2 )(2£ F ) 3 / 2 . The ana- 
lytic continuation by substitution ihw n — ► Hz — [i yields the 
asymptotic complex Greenfunction 



G(k, Z) M (-Hz)" 1 - £ k (-fe)- 2 + TTb(-hz) 



-5/2 



(A7) 



for \z\ — > oo. Eventually from (2.27) we obtain the asymp- 
totic spectral function A(k, e) « b 9(e) e~ 5 / 2 for e — > ±oo. 
Thus, the weight of the asymptotic power law of the spectral 
function is described by the constant factor b which is real, 
positive, and independent of k. 

A prior spectrum which meets all these requirements and 
which shows the correct wings is given by 



A (e) 



[e 2 + 7 



211/2 



2 [ e 2 +7 2]7/4 



(A8) 



The denominator represents a modified Lorentz spectrum with 
a non trivial exponent. In order to have a smooth function, we 
choose a large spectral width 7 = 20 £p. The specific form of 



the numerator and the exponent of the denominator guarantee 
the correct wings for e — > ±00 in leading order. It turns out 
that the prior spectrum Aq (e) can be chosen independent of 
k. 

In our implementation of the method we solve Eq. (A5) 
numerically by using Bryan's algorithm [77]. The rectangular 
matrix of the discretized constraint equations (2.17) is decom- 
posed by using a singular-value decomposition. Eventually 
we observe that only a small fraction of about 15-20 eigenval- 
ues provide essential contributions for the spectral function 
A(k,e). 

The maximum-entropy method must be applied for each 
value of the wave vector k in order to obtain the complete 
spectral function A(k, e). We find that the parameters of the 
method a, a, 7, and b can be chosen independent of k. We 
observe that the entropy (Al) together with the prior spectrum 
(A8) guarantees a positive spectral function A(k, e) > 0. Fi- 
nally, in our numerical calculations we observe that the di- 
mensionless difference (A3) has the same order ~ 10~ 5 over 
the whole range of Matsubara frequencies u n and for all k 
which is essential for the quality of our implementation of the 
maximum-entropy method. 



APPENDIX B: LIFETIME OF FERMIONIC EXCITATIONS 
AT ZERO TEMPERATURE 

In this Appendix we outline an analytical calculation of the 
lifetime of fermionic excitations at zero temperature, that is 
perturbative in deviations from the exactly soluble reduced 
BCS Hamiltonian (2.32). For arguments that indicate a break- 
down of well defined fermionic excitations in the opposite 
BEC limit see [78]. 

Quite generally a quasiparticle description of the BCS- 
BEC crossover problem requires that the low lying excitations 
above the exact ground state with energy Eq can be described 
by a non-interacting gas of quasiparticles 



H 



l k<7 



(Bl) 



The first term accounts for the Bogoliubo v- Anderson phonons 
with linear dispersion w q = c s q for momenta q that are small 
compared to the inverse healing length. The second term de- 
scribes the fermionic excitations which have a gapped spec- 
trum i?k- The crucial requirement that the lifetime of the 
quasiparticles by far exceeds their energy is trivially fulfilled 
for the bosonic excitations. In the weak coupling BCS-regime 
their lifetime is actually infinite up to an energy 2A, which is 
necessary for a decay into two fermionic quasiparticles. On 
the BEC-side they can decay through nonlinear corrections to 
the quantum hydrodynamic Lagrangian (3.2). Provided that 
the curvature parameter a introduced in Sec. Ill is negative, 
this leads to a width ~ q 5 by Beliaev damping, which is neg- 
ligible in the q — > limit. 

Regarding the fermionic quasiparticles, their lifetime turns 
out to be infinite near the dispersion minimum, despite the fact 
that their excitation energy becomes of the order of the Fermi 
energy in the experimentally relevant regime near a Feshbach 
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FIG. 7: Dominant contribution to the BCS-quasiparticle self energy 
at zero temperature. 



resonance. At unitarity, for instance, the zero temperature gap 
is A = 0.46c? f [21], in very good agreement with recent quan- 
tum Monte Carlo results [79]. The fermionic spectral func- 
tion should therefore exhibit a sharp peak near the dispersion 
minimum at zero temperature. Indeed, the relevant process 
that limits the lifetime close to the dispersion minimum is the 
emission of a Bogoliubov-Anderson phonon with momentum 
q. Due to energy- and momentum conservation, this process 
must obey the kinematic constraint 



c s |q| 



(B2) 



where k is the initial momentum of the fermionic excita- 
tion with dispersion E^ = /x + \J (ek — /i) 2 + A 2 and c s 
is the sound velocity. Equation (B2) implies, that the emis- 
sion of a phonon is impossible as long as the group velocity 
of the fermionic excitations is smaller than the sound veloc- 
ity \dEk/d\s\ < c s . This condition is always true for a small 
interval of momenta around the dispersion minimum, imply- 
ing that the lifetime of a fermionic excitation is infinite in this 
region. 

In the following we show briefly how the kinematic con- 
straint (B2) arises, if the residual interaction (2.33) between 
BCS quasiparticles is taken into account perturbatively. Our 
starting point is a reformulation of the BCS-BEC crossover 



Hamiltonian (2.1) in terms of BCS quasiparticle operators. 
The reduced BCS Hamiltonian (2.32) can be diagonalized ex- 
actly [43, 80] and takes the form 



H 



BCS 



TTlBCS 



E 

k.cr 



(B3) 



It has the form of the more general quasiparticle Hamilto- 
nian (Bl), but is actually valid at arbitrary momenta and 
energies. However, it misses completely the Bogoliubov- 
Anderson phonons. The quasiparticle operators are re- 
lated to the fermionic operators c kcr via the usual Bogoliubov 
transformation 



c k| 



'-kl = -"k Q kT + u k"-kl 



(B4) 
(B5) 



with the coefficients u| = [l + (e k — /i)/(_E k — fi))/2 and v 



[1 — (e k — fi)/(Ek — /u)]/2. The ground state is determined by 
the condition a kcr | GS) = 0. Before proceeding, we mention 
that the lifetime of the fermionic excitations is directly related 
to the lifetime of the BCS quasiparticles, since the fermionic 
Green's function G(k, u) can be expressed in terms of the 
BCS quasiparticle Green's function C? (k, oj) as 



G CT (k, lu) = u\ g a (k,uj) - vl <?-cr(-k, —u) 



(B6) 



using the Bogoliubov transformation defined in Eq. (B4) and 
(B5). Note that the second term in Eq. (B6) leads to the 
fermionic hole excitation spectrum with dispersion £? k ' (see 
Eq. (2.25)), even though the excitation energies of the BCS- 
quasiparticles are strictly positive. The residual interaction 
(2.33) describes the interaction between BCS quasiparticles 
and their coupling to the collective Bogoliubov-Anderson 
mode. Explicitely, Eq. (2.33) gives rise to three different types 
of quasiparticle interactions _ff les = H40+H31+H22 that have 
been discussed previously e.g. in nuclear physics [81] 



^ 40 = t7 E w k+Qt'k"k'Uk'+Q akTa-k-Qia-k'i a k'+QT + h.c. (B7) 

k,k',Q#0 

^ 31 = ^7 ( v k't'k' + QfkMk-Q - Uk'Wk' + QUkt'k-Q) aL a k-Q CT a -k'l a k' + QT + h - C - ( fi 8) 

k,k',Q#0,cr 

H22 = ^7 E [( M Q-kM k Uk'UQ-k' + fQ-kWkVk'VQ-k') a Q-kT a ki a k'i Q; Q-k'T 



k,k',Q#0 

+ (wkWk'-QWk'Wk-Q + «kVk'-Q^k'Wk-Q) a kT a -k'l a Q-k'l a k-QT 

+ Wk+QWk'+QVkVk' a k +QT a k'T a k'+QT a kT + Vk+QVk'+QWkWk' "k+Qi^k'i^k'+Qiaki 



(B9) 



corresponding to four-wave annihilation, quasiparticle decay the decay into three (or more) quasiparticles and the emission 
and quasiparticle scattering. The only processes that limit the of a Bogoliubov-Anderson phonon or a combination thereof, 
lifetime of a quasiparticle excitation at zero temperature are The decay into three quasiparticles via H 31 has a threshold 
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energy of 3 A and is forbidden in a rather broad range around 
the dispersion minimum. As discussed above, the emission 
of a Bogoliubov- Anderson phonon has a much less restrictive 
kinematic constraint and thus is the relevant lifetime-limiting 
process close to the dispersion minimum. In order to estimate 
this contribution, we need to know how the BCS quasiparti- 
cles couple to the collective Bogoliubov-Anderson mode. It is 
important to notice, that the phonons in the Hamiltonian (Bl) 
are not independent excitations, but are actually bound states 
of two BCS-quasiparticles. Indeed, as shown already by Gal- 



itskii [82], the vertex function T(q, u>) for the scattering of an 
up- and a down- BCS quasiparticle with total energy to and 
total momentum q has a pole at lo 2 = c? s q 2 corresponding to 
the Bogoliubov-Anderson phonon mode. Within a diagram- 
matic formulation, the leading order self-energy contribution 
to the BCS quasiparticle Green's function Q corresponding to 
the emission of a Bogoliubov-Anderson phonon is thus given 
by the diagram shown in Fig. 7. Explicitely, this gives rise to 
an imaginary part of the retarded self-energy given by 



Im^M = / V£&VS&> ^ba (?) 6{ U - £ k - q - cM) 

x Ree- ( r(k',w- J Bk- q - J B q -k')Ree- CT (k",w- J E; k _ q -S q _ k /0 + --. (BIO) 



where Zba(<7) denotes the quasiparticle weight of the 
Bogoliubov-Anderson mode and yt 1 - 3 ' is the bare vertex re- 
lated to H31. The dots indicate two more terms coming from 
the imaginary parts of the two BCS quasiparticle Green's 
functions in Eq. (BIO). However, these terms are not impor- 
tant close to the minimum of the dispersion curve, since they 
give rise to a kinematic constraint related to the decay into 
three quasiparticles which has a threshold energy of 3A. 

Assuming that the real part of the self-energy is small, we 
evaluate the self-energy at to = and extract from (BIO) the 
kinematic constraint (B2). Thus, in the weak coupling limit, 
the spectral function exhibits sharp peaks in an exponentially 
small interval 

\k — kp\ c, A 

' , ' < ^ (Bll) 

hp Lvp Bp 

around the minimum of the dispersion relation. Extrapolated 
to unitarity, the range where no broadening of the spectral 
function is expected is | k — | < O.lkp. In the BEC-regime, 



where the chemical potential is negative and the minimum of 
the dispersion relation is at k = 0, Eq. (B2) indicates that the 
spectral function is sharp for momenta k < mc s . 

Interestingly, the kinematic constraint (B2) that leads to an 
infinite lifetime of the fermionic excitations around the dis- 
persion minimum is implicitly also present in our Luttinger- 
Ward theory. Indeed, if Eq. (2.11) is reformulated in terms 
of mean-field Green's functions, the second term corresponds 
exactly to the self-energy contribution coming from the vir- 
tual emission of a Bogoliubov-Anderson phonon due to the 
phonon pole of the Vertex function 1\ Again, this process 
causes the constraint (B2). Nevertheless, our numerics show 
a finite lifetime at the dispersion minimum. Apart from the 
fact that a sharp feature in the spectral function can hardly be 
resolved numerically, we attribute the finite lifetime to the self 
consistent solution of the equations, since the replacement of 
bare with dressed Green's functions gives rise to diagrams that 
explicitly violate the Pauli principle. 
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